#opening csv
soils <- read_csv(here("data", "shift_soil_data.csv")) %>% # read in soil csv
mutate(date = lubridate::mdy(date), # change to date class
week = lubridate::week(date)) %>% # create a week column
select(soil_id:transect, meter_location, electro_cond_mS_per_cm, season, latitude, longitude)%>% #select certain column to work with
mutate(paired_flight = lubridate::ymd(paired_flight)) %>%
drop_na() %>% #drop any NA rows
st_as_sf(coords = c("longitude", "latitude")) %>%
st_set_crs(value = 4326) %>%
st_transform(crs = 32611)
# Open spatial Data
soils_sf <- read_sf(here("data", "soil_w_elevation.shp")) %>% # read in sf file
mutate(date = lubridate::ymd(date), #change date to date class
week = lubridate::week(date)) %>% #create week
select(electro_co, date, transect, meter_loca) %>% # select certain columns
filter(date != "2022-09-12") %>%
mutate(electro_co = as.numeric(electro_co)) %>% #change to numeric class
drop_na() #drop NAs
bbox <- read_sf(here("data", "bbox.shp")) #read in bounding box file
boundary <- read_sf(here("data", "boundary.shp")) %>%
st_transform(crs = 32611)
dem <- read_stars(here("data", "dem.tif")) %>% # read in dem
st_crop(y = st_bbox(bbox)) %>% #crop to bounding box
st_warp(cellsize = 4.8, crs = 32611)
mllw <- dem + 0.042
mllw_all <- c(mllw, mllw, mllw, mllw, mllw, mllw, mllw, mllw, mllw, mllw, mllw, mllw, mllw, along = 3) %>%
st_set_dimensions(3, values = lubridate::ymd(c("2022_02_24", "2022_02_28", "2022_03_08", "2022_03_16", "2022_03_22", "2022_04_05", "2022_04_12", "2022_04_20", "2022_04_29", "2022_05_03", "2022_05_11", "2022_05_17", "2022_05_29")), names = "date")
mllw_extract <- mllw %>% # call dem
st_extract(soils$geometry) # extract elevation to point layer
soils <- cbind(soils, mllw_extract$dem.tif) %>% #bind extracted values to the soil data
rename("elevation" = "mllw_extract.dem.tif") %>%
st_transform(crs = 32611)
#make descriptive stat table
soils %>% # to soils
group_by(transect) %>%# group by transect
summarize(min = min(electro_cond_mS_per_cm),#calculate mins
max = max(electro_cond_mS_per_cm),#calculate max
mean = mean(electro_cond_mS_per_cm),#calculate mean
dev = sd(electro_cond_mS_per_cm),#obtain standard dev
variance = var(electro_cond_mS_per_cm)) %>% # calculate variance
kableExtra::kable() %>% # create a table format
kableExtra::kable_classic(lightable_options = "striped")#theme the table
| transect | min | max | mean | dev | variance | geometry |
|---|---|---|---|---|---|---|
| COPR_1_EW | 0.42 | 1.34 | 0.6941176 | 0.2020538 | 0.0408257 | MULTIPOINT ((235669.2 38115… |
| COPR_2_EW | 6.32 | 10.53 | 8.2875000 | 1.2435749 | 1.5464786 | MULTIPOINT ((235855 3812030… |
| COPR_2_NS | 6.30 | 13.54 | 9.7581250 | 2.1347107 | 4.5569896 | MULTIPOINT ((235855 3812030… |
| NCOS_1_NS | 0.12 | 7.47 | 2.7527778 | 2.0644399 | 4.2619121 | MULTIPOINT ((235777.4 38125… |
| NCOS_2_EW | 0.74 | 8.62 | 2.9125714 | 1.8626205 | 3.4693550 | MULTIPOINT ((235671.2 38123… |
| NCOS_2_NS | 2.88 | 10.14 | 4.6538889 | 1.4925439 | 2.2276873 | MULTIPOINT ((235685.7 38123… |
# make descriptive stat dataframe
soil_stats <- soils %>% #call soils
group_by(transect) %>% #group by transects
summarize(mean = mean(electro_cond_mS_per_cm),# create dataframe with these stats
max = max(electro_cond_mS_per_cm),
min = min(electro_cond_mS_per_cm),
range = max - min)
ggplot(soils, aes(x = elevation, y = electro_cond_mS_per_cm, color = transect)) +#ggplot of soils/dem data
geom_point()+#geometry of the plot
scale_color_manual(values = calecopal::cal_palette(name = "superbloom3", n =6, type = "discrete"))+# colors to use
labs(x = "Elevation (m)",#labels
y = "Electrical conductivity (mS/cm)")+
theme(legend.position = "none")+#no legend theme
guides(color = guide_legend(title = "Transect ID"))+#changing legend title
ggtitle("Soil Electrical Conductivity Across Elevations")+# title of plot
theme(plot.title = element_text(color = "#5b4f41", size = 20, hjust = 0.5),# plot theming
plot.background = element_rect("white"),
panel.background = element_rect("#faf7f2"),
panel.grid = element_line(linetype= "longdash", color = "#f0ece1"),
axis.text = element_text(color = "#5b4f41", size = 14),
axis.title = element_text(color = "#5b4f41", size = 16),
strip.background = element_rect("#f8f8f8"),
strip.text = element_text(color = "#5b4f41"),
axis.line = element_line(color = "#5b4f41"),
legend.text = element_text(color = "#5b4f41"),
legend.title = element_text(color = "#5b4f41"))
strip_labels <- c("COPR_1_EW" = "COPR 1 EW",#create better formated labels for facet wrap strips
"COPR_2_EW" = "COPR 2 EW",
"COPR_2_NS" = "COPR 2 NS",
"NCOS_1_NS" = "NCOS 1 NS",
"NCOS_2_EW" = "NCOS 2 EW",
"NCOS_2_NS" = "NCOS 2 NS")
ggplot(soils, aes(x = elevation, y = electro_cond_mS_per_cm, color = transect)) +#start ggplot
geom_point()+# point geometry
scale_color_manual(values = calecopal::cal_palette(name = "superbloom3", n =6, type = "discrete"))+#colors
facet_wrap(~transect, labeller = as_labeller(strip_labels))+#subplots based on transect
theme(legend.position = "none")+# legend theme
labs(x = "Elevation (m)",#labels
y = "Electrical conductivity (mS/cm)")+
guides(color = guide_legend(title = "Transect ID"))+#legend title
ggtitle("Soil Electrical Conductivity Across Elevations")+#plot title
theme(plot.title = element_text(color = "#5b4f41", size = 16, hjust = 0.5),#plot theme
plot.background = element_rect("white"),
panel.background = element_rect("#faf7f2"),
panel.grid = element_line(linetype= "longdash", color = "#f0ece1"),
axis.text = element_text(color = "#5b4f41"),
axis.title = element_text(color = "#5b4f41"),
strip.background = element_rect("#f8f8f8"),
strip.text = element_text(color = "#5b4f41"),
axis.line = element_line(color = "#5b4f41"),
legend.text = element_text(color = "#5b4f41"),
legend.title = element_text(color = "#5b4f41"))
ndvi_1 <- read_stars(here("imagery", "NDVI", "2022_02_24_ndvi.tif")) %>%
setNames("ndvi")
ndvi_2 <- read_stars(here("imagery", "NDVI", "2022_02_28_ndvi.tif")) %>%
setNames("ndvi")
ndvi_3 <- read_stars(here("imagery", "NDVI", "2022_03_08_ndvi.tif")) %>%
setNames("ndvi")
ndvi_4 <- read_stars(here("imagery", "NDVI", "2022_03_16_ndvi.tif")) %>%
setNames("ndvi")
ndvi_5 <- read_stars(here("imagery", "NDVI", "2022_03_22_ndvi.tif")) %>%
setNames("ndvi")
ndvi_6 <- read_stars(here("imagery", "NDVI", "2022_04_05_ndvi.tif")) %>%
setNames("ndvi")
ndvi_7 <- read_stars(here("imagery", "NDVI", "2022_04_12_ndvi.tif")) %>%
setNames("ndvi")
ndvi_8 <- read_stars(here("imagery", "NDVI", "2022_04_20_ndvi.tif")) %>%
setNames("ndvi")
ndvi_9 <- read_stars(here("imagery", "NDVI", "2022_04_29_ndvi.tif")) %>%
setNames("ndvi")
ndvi_10 <- read_stars(here("imagery", "NDVI", "2022_05_03_ndvi.tif")) %>%
setNames("ndvi")
ndvi_11 <- read_stars(here("imagery", "NDVI", "2022_05_11_ndvi.tif")) %>%
setNames("ndvi")
#ndvi_12 <- read_stars(here("imagery", "NDVI", "2022_05_12_ndvi.tif")) %>%
# setNames("ndvi")
ndvi_13 <- read_stars(here("imagery", "NDVI", "2022_05_17_ndvi.tif")) %>%
setNames("ndvi")
ndvi_14 <- read_stars(here("imagery", "NDVI", "2022_05_29_ndvi.tif")) %>%
setNames("ndvi")
ndvi_all <- c(ndvi_1, ndvi_2, ndvi_3, ndvi_4, ndvi_5, ndvi_6, ndvi_7, ndvi_8, ndvi_9, ndvi_10, ndvi_11, ndvi_13, ndvi_14, along = 3) %>%
st_set_dimensions(3, values = lubridate::ymd(c("2022_02_24", "2022_02_28", "2022_03_08", "2022_03_16", "2022_03_22", "2022_04_05", "2022_04_12", "2022_04_20", "2022_04_29", "2022_05_03", "2022_05_11", "2022_05_17", "2022_05_29")), names = "date")
vssi_1 <- read_stars(here("imagery", "VSSI", "02_24_vssi.tif")) %>%
setNames("vssi")
vssi_2 <- read_stars(here("imagery", "VSSI", "02_28_vssi.tif")) %>%
setNames("vssi")
vssi_3 <- read_stars(here("imagery", "VSSI", "03_08_vssi.tif")) %>%
setNames("vssi")
vssi_4 <- read_stars(here("imagery", "VSSI", "03_16_vssi.tif")) %>%
setNames("vssi")
vssi_5 <- read_stars(here("imagery", "VSSI", "03_22_vssi.tif")) %>%
setNames("vssi")
vssi_6 <- read_stars(here("imagery", "VSSI", "04_05_vssi.tif")) %>%
setNames("vssi")
vssi_7 <- read_stars(here("imagery", "VSSI", "04_12_vssi.tif")) %>%
setNames("vssi")
vssi_8 <- read_stars(here("imagery", "VSSI", "04_20_vssi.tif")) %>%
setNames("vssi")
vssi_9 <- read_stars(here("imagery", "VSSI", "04_29_vssi.tif")) %>%
setNames("vssi")
vssi_10 <- read_stars(here("imagery", "VSSI", "05_03_vssi.tif")) %>%
setNames("vssi")
vssi_11 <- read_stars(here("imagery", "VSSI", "05_11_vssi.tif")) %>%
setNames("vssi")
#vssi_12 <- read_stars(here("imagery", "VSSI", "05_12_vssi.tif")) %>%
# setNames("vssi")
vssi_13 <- read_stars(here("imagery", "VSSI", "05_17_vssi.tif")) %>%
setNames("vssi")
vssi_14 <- read_stars(here("imagery", "VSSI", "05_29_vssi")) %>%
setNames("vssi")
vssi_all <- c(vssi_1, vssi_2, vssi_3, vssi_4, vssi_5, vssi_6, vssi_7, vssi_8, vssi_9, vssi_10, vssi_11, vssi_13, vssi_14, along = 3) %>%
st_set_dimensions(3, values = lubridate::ymd(c("2022_02_24", "2022_02_28", "2022_03_08", "2022_03_16", "2022_03_22", "2022_04_05", "2022_04_12", "2022_04_20", "2022_04_29", "2022_05_03", "2022_05_11", "2022_05_17", "2022_05_29")), names = "date")
mari_1 <- read_stars(here("imagery", "mARI", "02_24_mari_avg.tif")) %>%
setNames("mari")
mari_2 <- read_stars(here("imagery", "mARI", "02_28_mari_avg.tif")) %>%
setNames("mari")
mari_3 <- read_stars(here("imagery", "mARI", "03_08_mari_avg.tif")) %>%
setNames("mari")
mari_4 <- read_stars(here("imagery", "mARI", "03_16_mari_avg.tif")) %>%
setNames("mari")
mari_5 <- read_stars(here("imagery", "mARI", "03_22_mari_avg.tif")) %>%
setNames("mari")
mari_6 <- read_stars(here("imagery", "mARI", "04_05_mari_avg.tif")) %>%
setNames("mari")
mari_7 <- read_stars(here("imagery", "mARI", "04_12_mari_avg.tif")) %>%
setNames("mari")
mari_8 <- read_stars(here("imagery", "mARI", "04_20_mari_avg.tif")) %>%
setNames("mari")
mari_9 <- read_stars(here("imagery", "mARI", "04_29_mari_avg.tif")) %>%
setNames("mari")
mari_10 <- read_stars(here("imagery", "mARI", "05_03_mari_avg.tif")) %>%
setNames("mari")
mari_11 <- read_stars(here("imagery", "mARI", "05_11_mari_avg.tif")) %>%
setNames("mari")
#mari_12 <- read_stars(here("imagery", "mARI", "05_12_mari_avg.tif")) %>%
# setNames("mari")
mari_13 <- read_stars(here("imagery", "mARI", "05_17_mari_avg.tif")) %>%
setNames("mari")
mari_14 <- read_stars(here("imagery", "mARI", "05_29_mari_avg.tif")) %>%
setNames("mari")
mari_all <- c(mari_1, mari_2, mari_3, mari_4, mari_5, mari_6, mari_7, mari_8, mari_9, mari_10, mari_11, mari_13, mari_14, along = 3) %>%
st_set_dimensions(3, values = lubridate::ymd(c("2022_02_24", "2022_02_28", "2022_03_08", "2022_03_16", "2022_03_22", "2022_04_05", "2022_04_12", "2022_04_20", "2022_04_29", "2022_05_03", "2022_05_11", "2022_05_17", "2022_05_29")), names = "date")
crsi_1 <- read_stars(here("imagery", "CRSI", "02_24_crsi.tif")) %>%
setNames("crsi")
crsi_2 <- read_stars(here("imagery", "CRSI", "02_28_crsi.tif")) %>%
setNames("crsi")
crsi_3 <- read_stars(here("imagery", "CRSI", "03_08_crsi.tif")) %>%
setNames("crsi")
crsi_4 <- read_stars(here("imagery", "CRSI", "03_16_crsi.tif")) %>%
setNames("crsi")
crsi_5 <- read_stars(here("imagery", "CRSI", "03_22_crsi.tif")) %>%
setNames("crsi")
crsi_6 <- read_stars(here("imagery", "CRSI", "04_05_crsi.tif")) %>%
setNames("crsi")
crsi_7 <- read_stars(here("imagery", "CRSI", "04_12_crsi.tif")) %>%
setNames("crsi")
crsi_8 <- read_stars(here("imagery", "CRSI", "04_20_crsi.tif")) %>%
setNames("crsi")
crsi_9 <- read_stars(here("imagery", "CRSI", "04_29_crsi.tif")) %>%
setNames("crsi")
crsi_10 <- read_stars(here("imagery", "CRSI", "05_03_crsi.tif")) %>%
setNames("crsi")
crsi_11 <- read_stars(here("imagery", "CRSI", "05_11_crsi.tif")) %>%
setNames("crsi")
#crsi_12 <- read_stars(here("imagery", "CRSI", "05_12_crsi.tif")) %>%
# setNames("crsi")
crsi_13 <- read_stars(here("imagery", "CRSI", "05_17_crsi.tif")) %>%
setNames("crsi")
crsi_14 <- read_stars(here("imagery", "CRSI", "05_29_crsi.tif")) %>%
setNames("crsi")
crsi_all <- c(crsi_1, crsi_2, crsi_3, crsi_4, crsi_5, crsi_6, crsi_7, crsi_8, crsi_9, crsi_10, crsi_11, crsi_13, crsi_14, along = 3) %>%
st_set_dimensions(3, values = lubridate::ymd(c("2022_02_24", "2022_02_28", "2022_03_08", "2022_03_16", "2022_03_22", "2022_04_05", "2022_04_12", "2022_04_20", "2022_04_29", "2022_05_03", "2022_05_11", "2022_05_17", "2022_05_29")), names = "date")
nir_1 <- read_stars(here("imagery", "subset_images", "rfl_2022_02_24.tif")) %>%
st_set_dimensions(3, values = seq(1:425)) %>%
filter(band %in% c(76:126)) %>%
split(3) %>%
setNames(c(76:126))
nir_2 <- read_stars(here("imagery", "subset_images", "rfl_2022_02_28.tif")) %>%
st_set_dimensions(3, values = seq(1:425)) %>%
filter(band %in% c(76:126)) %>%
split(3) %>%
setNames(c(76:126))
nir_3 <- read_stars(here("imagery", "subset_images", "rfl_2022_03_08.tif")) %>%
st_set_dimensions(3, values = seq(1:425)) %>%
filter(band %in% c(76:126)) %>%
split(3) %>%
setNames(c(76:126))
nir_4 <- read_stars(here("imagery", "subset_images", "rfl_2022_03_16.tif")) %>%
st_set_dimensions(3, values = seq(1:425)) %>%
filter(band %in% c(76:126)) %>%
split(3) %>%
setNames(c(76:126))
nir_5 <- read_stars(here("imagery", "subset_images", "rfl_2022_03_22.tif")) %>%
st_set_dimensions(3, values = seq(1:425)) %>%
filter(band %in% c(76:126)) %>%
split(3) %>%
setNames(c(76:126))
nir_6 <- read_stars(here("imagery", "subset_images", "rfl_2022_04_05.tif")) %>%
st_set_dimensions(3, values = seq(1:425)) %>%
filter(band %in% c(76:126)) %>%
split(3) %>%
setNames(c(76:126))
nir_7 <- read_stars(here("imagery", "subset_images", "rfl_2022_04_12.tif")) %>%
st_set_dimensions(3, values = seq(1:425)) %>%
filter(band %in% c(76:126)) %>%
split(3) %>%
setNames(c(76:126))
nir_8 <- read_stars(here("imagery", "subset_images", "rfl_2022_04_20.tif")) %>%
st_set_dimensions(3, values = seq(1:425)) %>%
filter(band %in% c(76:126)) %>%
split(3) %>%
setNames(c(76:126))
nir_9 <- read_stars(here("imagery", "subset_images", "rfl_2022_04_29.tif")) %>%
st_set_dimensions(3, values = seq(1:425)) %>%
filter(band %in% c(76:126)) %>%
split(3) %>%
setNames(c(76:126))
nir_10 <- read_stars(here("imagery", "subset_images", "rfl_2022_05_03.tif")) %>%
st_set_dimensions(3, values = seq(1:425)) %>%
filter(band %in% c(76:126)) %>%
split(3) %>%
setNames(c(76:126))
nir_11 <- read_stars(here("imagery", "subset_images", "rfl_2022_05_11.tif")) %>%
st_set_dimensions(3, values = seq(1:425)) %>%
filter(band %in% c(76:126)) %>%
split(3) %>%
setNames(c(76:126))
#nir_12 <- read_stars(here("imagery", "subset_images", "rfl_2022_05_12.tif")) %>%
# st_set_dimensions(3, values = seq(1:425)) %>%
# filter(band %in% c(76:126)) %>%
# split(3) %>%
# setNames(c(76:126))
nir_13 <- read_stars(here("imagery", "subset_images", "rfl_2022_05_17.tif")) %>%
st_set_dimensions(3, values = seq(1:425)) %>%
filter(band %in% c(76:126)) %>%
split(3) %>%
setNames(c(76:126))
nir_14 <- read_stars(here("imagery", "subset_images", "rfl_2022_05_29.tif")) %>%
st_set_dimensions(3, values = seq(1:425)) %>%
filter(band %in% c(76:126)) %>%
split(3) %>%
setNames(c(76:126))
nir_all <- c(nir_1, nir_2, nir_3, nir_4, nir_5, nir_6, nir_7, nir_8, nir_9, nir_10, nir_11, nir_13, nir_14, along = 3) %>%
st_set_dimensions(3, values = lubridate::ymd(c("2022_02_24", "2022_02_28", "2022_03_08", "2022_03_16", "2022_03_22", "2022_04_05", "2022_04_12", "2022_04_20", "2022_04_29", "2022_05_03", "2022_05_11", "2022_05_17", "2022_05_29")), names = "date") %>%
st_warp(dest = crsi_all)
mllw_all <- mllw_all %>%
st_warp(dest = crsi_all) %>%
setNames("elevation")
soils_extract <- data.frame()
for (i in unique(lubridate::week(soils$paired_flight))){
ndvi_week <- ndvi_all %>%
filter(lubridate::week(date) == i)
mari_week <- mari_all %>%
filter(lubridate::week(date) == i)
vssi_week <- vssi_all %>%
filter(lubridate::week(date) == i)
crsi_week <- crsi_all %>%
filter(lubridate::week(date) == i)
nir_week <- nir_all %>%
filter(lubridate::week(date) == i) %>%
st_as_sf
soil_week <- soils %>%
mutate(week = lubridate::week(paired_flight)) %>%
filter(week == i)
ndvi_extract <- ndvi_week %>%
st_extract(soil_week)
mari_extract <- mari_week %>%
st_extract(soil_week)
vssi_extract <- vssi_week %>%
st_extract(soil_week)
crsi_extract <- crsi_week %>%
st_extract(soil_week)
nir_extract <- soil_week %>%
st_join(nir_week) %>%
select("76":"126") %>%
st_drop_geometry()
soils_week <- bind_cols(soil_week, ndvi_extract$ndvi, mari_extract$mari, vssi_extract$vssi, crsi_extract$crsi, nir_extract) %>%
rename("ndvi" = "...11",
"mari" = "...12",
"vssi" = "...13",
"crsi" = "...14")
soils_extract <- rbind(soils_extract, soils_week)
}
# Read in the data
nir_pcd <- soils_extract %>% # call object
as.data.frame() %>%
select(electro_cond_mS_per_cm,
elevation,
ndvi:crsi,
`76`:`126`) %>% # select the variables we want to explore
drop_na() # drops columns with NAs as PCA only works with numeric values
# Performing the PCA
nir_pca <- nir_pcd %>% # pipe in the data
prcomp(scale. = TRUE) # rescaling of data is performed
summary(nir_pca)
## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 7.1759 1.75884 1.06655 0.77939 0.64010 0.36867 0.22789
## Proportion of Variance 0.9034 0.05427 0.01996 0.01066 0.00719 0.00238 0.00091
## Cumulative Proportion 0.9034 0.95767 0.97763 0.98829 0.99548 0.99786 0.99877
## PC8 PC9 PC10 PC11 PC12 PC13 PC14
## Standard deviation 0.17019 0.08872 0.07819 0.06710 0.05182 0.04731 0.04508
## Proportion of Variance 0.00051 0.00014 0.00011 0.00008 0.00005 0.00004 0.00004
## Cumulative Proportion 0.99928 0.99942 0.99952 0.99960 0.99965 0.99969 0.99973
## PC15 PC16 PC17 PC18 PC19 PC20 PC21
## Standard deviation 0.03938 0.03687 0.03189 0.03122 0.02883 0.02796 0.02572
## Proportion of Variance 0.00003 0.00002 0.00002 0.00002 0.00001 0.00001 0.00001
## Cumulative Proportion 0.99975 0.99978 0.99979 0.99981 0.99983 0.99984 0.99985
## PC22 PC23 PC24 PC25 PC26 PC27 PC28
## Standard deviation 0.02501 0.02379 0.02299 0.02248 0.02156 0.02062 0.02012
## Proportion of Variance 0.00001 0.00001 0.00001 0.00001 0.00001 0.00001 0.00001
## Cumulative Proportion 0.99986 0.99987 0.99988 0.99989 0.99990 0.99991 0.99991
## PC29 PC30 PC31 PC32 PC33 PC34 PC35
## Standard deviation 0.01975 0.01867 0.01835 0.01813 0.01700 0.0165 0.01604
## Proportion of Variance 0.00001 0.00001 0.00001 0.00001 0.00001 0.0000 0.00000
## Cumulative Proportion 0.99992 0.99993 0.99993 0.99994 0.99994 1.0000 0.99995
## PC36 PC37 PC38 PC39 PC40 PC41 PC42
## Standard deviation 0.01574 0.01485 0.01443 0.01417 0.01311 0.01259 0.0125
## Proportion of Variance 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.0000
## Cumulative Proportion 0.99996 0.99996 0.99996 0.99997 0.99997 0.99997 1.0000
## PC43 PC44 PC45 PC46 PC47 PC48 PC49
## Standard deviation 0.01182 0.01164 0.01121 0.01095 0.01063 0.01046 0.009629
## Proportion of Variance 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.000000
## Cumulative Proportion 0.99998 0.99998 0.99998 0.99999 0.99999 0.99999 0.999990
## PC50 PC51 PC52 PC53 PC54 PC55
## Standard deviation 0.009267 0.009191 0.008831 0.00823 0.007837 0.006979
## Proportion of Variance 0.000000 0.000000 0.000000 0.00000 0.000000 0.000000
## Cumulative Proportion 0.999990 0.999990 1.000000 1.00000 1.000000 1.000000
## PC56 PC57
## Standard deviation 0.006765 0.006138
## Proportion of Variance 0.000000 0.000000
## Cumulative Proportion 1.000000 1.000000
#Creating the biplot via autoplot()
autoplot(nir_pca, # pca data
x = 1,
y = 2,
data = nir_pcd, # original data
loadings = TRUE, # show loadings
loadings.label = TRUE, # label loadings
loadings.colour = "black", # loading color
loadings.label.colour = "black", # loading label color
loadings.label.vjust = -.5) + # vertical justification
ggtitle("Principle Component 1 and 2 Biplot") + # title
labs(x = "Principle Component 1",
y= "Principle Component 2")+ # axis labels
theme(plot.title = element_text(color = "#5b4f41", size = 16),
plot.background = element_rect("white"),
panel.background = element_rect("#faf7f2"),
panel.grid = element_line(linetype= "longdash",
color = "#f0ece1"),
axis.text = element_text(color = "#5b4f41"),
axis.title = element_text(color = "#5b4f41"),
strip.background = element_rect("white"),
axis.line = element_line(color = "#5b4f41")) # change themes
autoplot(nir_pca, # pca data
x = 1,
y = 3,
data = nir_pcd, # original data
loadings = TRUE, # show loadings
loadings.label = TRUE, # label loadings
loadings.colour = "black", # loading color
loadings.label.colour = "black", # loading label color
loadings.label.vjust = -.5) + # vertical justification
ggtitle("Principle Component 1 and 3 Biplot") + # title
labs(x = "Principle Component 1",
y= "Principle Component 3")+ # axis labels
theme(plot.title = element_text(color = "#5b4f41", size = 16),
plot.background = element_rect("white"),
panel.background = element_rect("#faf7f2"),
panel.grid = element_line(linetype= "longdash",
color = "#f0ece1"),
axis.text = element_text(color = "#5b4f41"),
axis.title = element_text(color = "#5b4f41"),
strip.background = element_rect("white"),
axis.line = element_line(color = "#5b4f41")) # change themes
autoplot(nir_pca, # pca data
x = 1,
y = 4,
data = nir_pcd, # original data
loadings = TRUE, # show loadings
loadings.label = TRUE, # label loadings
loadings.colour = "black", # loading color
loadings.label.colour = "black", # loading label color
loadings.label.vjust = -.5) + # vertical justification
ggtitle("Principle Component 1 and 4 Biplot") + # title
labs(x = "Principle Component 1",
y= "Principle Component 4")+ # axis labels
theme(plot.title = element_text(color = "#5b4f41", size = 16),
plot.background = element_rect("white"),
panel.background = element_rect("#faf7f2"),
panel.grid = element_line(linetype= "longdash",
color = "#f0ece1"),
axis.text = element_text(color = "#5b4f41"),
axis.title = element_text(color = "#5b4f41"),
strip.background = element_rect("white"),
axis.line = element_line(color = "#5b4f41")) # change themes
autoplot(nir_pca, # pca data
x = 2,
y = 3,
data = nir_pcd, # original data
loadings = TRUE, # show loadings
loadings.label = TRUE, # label loadings
loadings.colour = "black", # loading color
loadings.label.colour = "black", # loading label color
loadings.label.vjust = -.5) + # vertical justification
ggtitle("Principle Component 2 and 3 Biplot") + # title
labs(x = "Principle Component 2",
y= "Principle Component 3")+ # axis labels
theme(plot.title = element_text(color = "#5b4f41", size = 16),
plot.background = element_rect("white"),
panel.background = element_rect("#faf7f2"),
panel.grid = element_line(linetype= "longdash",
color = "#f0ece1"),
axis.text = element_text(color = "#5b4f41"),
axis.title = element_text(color = "#5b4f41"),
strip.background = element_rect("white"),
axis.line = element_line(color = "#5b4f41")) # change themes
autoplot(nir_pca, # pca data
x = 2,
y = 4,
data = nir_pcd, # original data
loadings = TRUE, # show loadings
loadings.label = TRUE, # label loadings
loadings.colour = "black", # loading color
loadings.label.colour = "black", # loading label color
loadings.label.vjust = -.5) + # vertical justification
ggtitle("Principle Component 2 and 4 Biplot") + # title
labs(x = "Principle Component 2",
y= "Principle Component 4")+ # axis labels
theme(plot.title = element_text(color = "#5b4f41", size = 16),
plot.background = element_rect("white"),
panel.background = element_rect("#faf7f2"),
panel.grid = element_line(linetype= "longdash",
color = "#f0ece1"),
axis.text = element_text(color = "#5b4f41"),
axis.title = element_text(color = "#5b4f41"),
strip.background = element_rect("white"),
axis.line = element_line(color = "#5b4f41")) # change themes
autoplot(nir_pca, # pca data
x = 3,
y = 4,
data = nir_pcd, # original data
loadings = TRUE, # show loadings
loadings.label = TRUE, # label loadings
loadings.colour = "black", # loading color
loadings.label.colour = "black", # loading label color
loadings.label.vjust = -.5) + # vertical justification
ggtitle("Principle Component 3 and 4 Biplot") + # title
labs(x = "Principle Component 3",
y= "Principle Component 4")+ # axis labels
theme(plot.title = element_text(color = "#5b4f41", size = 16),
plot.background = element_rect("white"),
panel.background = element_rect("#faf7f2"),
panel.grid = element_line(linetype= "longdash",
color = "#f0ece1"),
axis.text = element_text(color = "#5b4f41"),
axis.title = element_text(color = "#5b4f41"),
strip.background = element_rect("white"),
axis.line = element_line(color = "#5b4f41")) # change themes
# Creating screeplot
sd_vec <- nir_pca$sdev # standard deviation
var_vec <- sd_vec^2 # variance
pc_names <- colnames(nir_pca$rotation) # names of PCs
pct_expl_df <- data.frame(v = var_vec, # new data frame for screeplot; variance
pct_v = var_vec / sum(var_vec), # percent of variance
pc = fct_inorder(pc_names)) %>% # orders rows
mutate(pct_label = paste0(round(pct_v * 100, 1), "%")) # adds character %
ggplot(pct_expl_df, aes(x = pc, y = v))+ # graphs via PC and variance
geom_col()+
geom_text(aes(label = pct_label), vjust = 0, nudge_y = 0.008)+
labs(x = "Principle Components",
y = "Variances") +
ggtitle("PCA Screeplot")+
theme(plot.title = element_text(hjust = .5))
fit_all <- train(electro_cond_mS_per_cm ~ elevation + ndvi + mari + vssi + crsi,
data = soils_extract,
method = "rf",
trControl = trainControl(method = "cv", number = 10),
importance = TRUE)
# using caret
fit_all_imp <- as.data.frame(fit_all$finalModel$importance)
fit_all_imp_scaled <- predict(preProcess(fit_all_imp, method = c("range")), fit_all_imp) %>%
rownames_to_column(var = "variable")
ggplot(fit_all_imp_scaled) +
geom_point(aes(x = IncNodePurity, y = variable), color = "#E69512")+
geom_segment(aes(y= variable, x = 0, yend = variable, xend = IncNodePurity), color = "#E69512", size = 1.2)+
geom_point(aes(x = fit_all_imp_scaled$`%IncMSE`, y = variable), color ="#3A5D3D")+
geom_segment(aes(y= variable, x = 0, yend = variable, xend = fit_all_imp_scaled$`%IncMSE`), color = "#3A5D3D", linetype = "dashed")+
labs(y = "Variables",#labels
x = "Importance (Scaled)")+
ggtitle("Random Forest Variable Importance")+#plot title
theme(plot.title = element_text(color = "#5b4f41", size = 16, hjust = 0.5),#plot themes
plot.background = element_rect("white"),
panel.background = element_rect("#faf7f2"),
panel.grid = element_line(linetype= "longdash", color = "#f0ece1"),
axis.text = element_text(color = "#5b4f41"),
axis.title = element_text(color = "#5b4f41"),
strip.background = element_rect("#f8f8f8"),
strip.text = element_text(color = "#5b4f41"),
axis.line = element_line(color = "#5b4f41"))
cor(fit_all$finalModel$y, fit_all$finalModel$predicted)
## [1] 0.8421477
ggplot()+#ggplot object
geom_abline(slope = 1, intercept = 0, linetype = "dotted")+
geom_point(aes(y = fit_all$finalModel$y, x = fit_all$finalModel$predicted), color = "#E69512")+#data to make points from
geom_smooth(method = "lm", aes(y= fit_all$finalModel$y, x = fit_all$finalModel$predicted), color = "#3A5D3D")+#smooth line data
scale_x_continuous(expand = c(0, 0), limits = c(0, 14)) +
scale_y_continuous(expand = c(0, 0), limits = c(0, 14))+
labs(y = "Actual Electrical Conductivity (mS/cm)",#labels
x = "Predicted Electrical conductivity (mS/cm)")+
ggtitle("Elevation and Indices (caret method)")+#plot title
theme(plot.title = element_text(color = "#5b4f41", size = 16, hjust = 0.5),#plot themes
plot.background = element_rect("white"),
panel.background = element_rect("#faf7f2"),
panel.grid = element_line(linetype= "longdash", color = "#f0ece1"),
axis.text = element_text(color = "#5b4f41"),
axis.title = element_text(color = "#5b4f41"),
strip.background = element_rect("#f8f8f8"),
strip.text = element_text(color = "#5b4f41"),
axis.line = element_line(color = "#5b4f41"))
#ggsave(filename = "all_plot_caret.png", plot = last_plot(), width = 12, height = 8, device = "png", dpi = 500) #write an image of the plot
all_corr <- cor(y= fit_all$finalModel$y, x = fit_all$finalModel$predicted)
all_corr
## [1] 0.8421477
fit_si <- train(electro_cond_mS_per_cm ~ ndvi + mari + vssi + crsi + `76` + `77` + `78` + `79` + `80` + `81` + `82` + `83` + `84` + `85` + `86` + `87` + `88` + `89` + `90` + `91` + `92` + `93` + `94` + `95` + `96` + `97` + `98` + `99` + `100` + `101` + `102` + `103` + `104` + `105` + `106` + `107` + `108` + `109` + `110` + `111` + `112` + `113` + `114` + `115` + `116` + `117` + `118` + `119` + `120` + `121` + `122` + `123` + `124` + `125` + `126`,
data = soils_extract,
method = "rf",
trControl = trainControl(method = "cv", number = 10),
importance = TRUE)
fit_si_imp <- as.data.frame(fit_si$finalModel$importance)
fit_si_imp_scaled <- predict(preProcess(fit_si_imp, method = c("range")), fit_si_imp) %>%
rownames_to_column(var = "variable")
ggplot(fit_si_imp_scaled) +
geom_point(aes(x = IncNodePurity, y = variable), color = "#E69512")+
geom_segment(aes(y= variable, x = 0, yend = variable, xend = IncNodePurity), color = "#E69512", size = 1.2)+
geom_point(aes(x = fit_si_imp_scaled$`%IncMSE`, y = variable), color ="#3A5D3D")+
geom_segment(aes(y= variable, x = 0, yend = variable, xend = fit_si_imp_scaled$`%IncMSE`), color = "#3A5D3D", linetype = "dashed")+
labs(y = "Variables",#labels
x = "Importance (Scaled)")+
ggtitle("Random Forest Variable Importance")+#plot title
theme(plot.title = element_text(color = "#5b4f41", size = 16, hjust = 0.5),#plot themes
plot.background = element_rect("white"),
panel.background = element_rect("#faf7f2"),
panel.grid = element_line(linetype= "longdash", color = "#f0ece1"),
axis.text = element_text(color = "#5b4f41"),
axis.title = element_text(color = "#5b4f41"),
strip.background = element_rect("#f8f8f8"),
strip.text = element_text(color = "#5b4f41"),
axis.line = element_line(color = "#5b4f41"))
cor(y = fit_si$finalModel$y, x= fit_si$finalModel$predicted)
## [1] 0.5060143
ggplot()+#ggplot object
geom_abline(slope = 1, intercept = 0, linetype = "dotted")+
geom_point(aes(y= fit_si$finalModel$y, x = fit_si$finalModel$predicted), color = "#E69512")+#data to make points from
geom_smooth(method = "lm", aes(y= fit_si$finalModel$y, x = fit_si$finalModel$predicted), color = "#3A5D3D")+#smooth line data
scale_x_continuous(expand = c(0, 0), limits = c(0, 14)) +
scale_y_continuous(expand = c(0, 0), limits = c(0, 14))+
labs(y = "Actual Electrical Conductivity (mS/cm)",#labels
x = "Predicted Electrical conductivity (mS/cm)")+
ggtitle("Spectral Indices")+#plot title
theme(plot.title = element_text(color = "#5b4f41", size = 16, hjust = 0.5),#plot themes
plot.background = element_rect("white"),
panel.background = element_rect("#faf7f2"),
panel.grid = element_line(linetype= "longdash", color = "#f0ece1"),
axis.text = element_text(color = "#5b4f41"),
axis.title = element_text(color = "#5b4f41"),
strip.background = element_rect("#f8f8f8"),
strip.text = element_text(color = "#5b4f41"),
axis.line = element_line(color = "#5b4f41"))
#ggsave(filename = "si_plot.png", plot = last_plot(), width = 12, height = 8, device = "png", dpi = 500) #write an image of the plot
fit_ele <- train(electro_cond_mS_per_cm ~ elevation, data = soils, method = "rf", trControl = trainControl(method = "cv", number = 10), importance = T)#random forest fitting of the data with stated equation
elevation_cor <- cor(y = fit_ele$finalModel$y, x = fit_ele$finalModel$predicted)#calculate correlation coefficient
elevation_cor#call/print correlation value
## [1] 0.9208296
ggplot()+#ggplot object
geom_abline(slope = 1, intercept = 0, linetype = "dotted")+
geom_point(aes(y = fit_ele$finalModel$y, x = fit_ele$finalModel$predicted), color = "#E69512")+#data to make points from
geom_smooth(method = "lm", aes(y = fit_ele$finalModel$y, x = fit_ele$finalModel$predicted), color = "#3A5D3D")+#smooth line data
scale_x_continuous(expand = c(0, 0), limits = c(0, 14)) +
scale_y_continuous(expand = c(0, 0), limits = c(0, 14))+
labs(y = "Actual Electrical Conductivity (mS/cm)",#labels
x = "Predicted Electrical conductivity (mS/cm)")+
ggtitle("Elevation")+#plot title
theme(plot.title = element_text(color = "#5b4f41", size = 16, hjust = 0.5),#plot themes
plot.background = element_rect("white"),
panel.background = element_rect("#faf7f2"),
panel.grid = element_line(linetype= "longdash", color = "#f0ece1"),
axis.text = element_text(color = "#5b4f41"),
axis.title = element_text(color = "#5b4f41"),
strip.background = element_rect("#f8f8f8"),
strip.text = element_text(color = "#5b4f41"),
axis.line = element_line(color = "#5b4f41"))
#ggsave(filename = "ele_plot.png", plot = last_plot(), width = 12, height = 8, device = "png", dpi = 500) #write an image of the plot
raster_all <- c(crsi_all, ndvi_all, vssi_all, mari_all, mllw_all)
all_pred <- predict(raster_all, fit_all, drop_dimensions = F)
all_pred_feb <- all_pred %>%
filter(date == "2022-02-24") %>%
st_crop(y = boundary, crop = TRUE, as_points = FALSE)
ggplot()+
geom_stars(data = all_pred_feb, aes(x, y, fill = prediction)) +
scale_fill_gradientn(colors = rev(calecopal::cal_palette("desert", 2, "continuous")))+
labs(y = "Northing",#labels
x = "Easting",
fill = "Electrical
Conductivity
(mS/cm)")+
ggtitle("Soil Salinity")+#plot title
theme(plot.title = element_text(color = "#5b4f41", size = 16, hjust = 0.5),#plot themes
plot.background = element_rect("white"),
panel.background = element_rect("#faf7f2"),
panel.grid = element_line(linetype= "longdash", color = "#f0ece1"),
axis.text = element_text(color = "#5b4f41"),
axis.title = element_text(color = "#5b4f41"),
strip.background = element_rect("#f8f8f8"),
strip.text = element_text(color = "#5b4f41"),
axis.line = element_line(color = "#5b4f41"),
legend.title = element_text(color = "#5b4f41"),
legend.text = element_text(color = "#5b4f41"))
# Interactive Map
mapview::mapview(all_pred_feb)
fit_all_pls <- train(electro_cond_mS_per_cm ~ elevation + ndvi + mari + vssi + crsi,
data = soils_extract,
method = "pls",
trControl = trainControl(method = "cv", number = 10),
importance = TRUE)
# using caret
predicted_values <- as.data.frame(fit_all_pls$finalModel$fitted.values) %>%
select(`.outcome.3 comps`)
training_data <- as.data.frame(fit_all_pls$trainingData$.outcome)
train_pred <- c(predicted_values, training_data) %>%
as.data.frame()
cor(training_data, predicted_values)
## .outcome.3 comps
## fit_all_pls$trainingData$.outcome 0.6955281
ggplot()+#ggplot object
geom_abline(slope = 1, intercept = 0, linetype = "dotted")+
geom_point(data = train_pred, aes(y = fit_all_pls.trainingData..outcome, x = `.outcome.3.comps`), color = "#E69512")+#data to make points from
geom_smooth(method = "lm", data = train_pred, aes(y = fit_all_pls.trainingData..outcome, x = `.outcome.3.comps`), color = "#3A5D3D")+#smooth line data
labs(y = "Actual Electrical Conductivity (mS/cm)",#labels
x = "Predicted Electrical conductivity (mS/cm)")+
ggtitle("Elevation and Indices (PLSR)")+#plot title
theme(plot.title = element_text(color = "#5b4f41", size = 16, hjust = 0.5),#plot themes
plot.background = element_rect("white"),
panel.background = element_rect("#faf7f2"),
panel.grid = element_line(linetype= "longdash", color = "#f0ece1"),
axis.text = element_text(color = "#5b4f41"),
axis.title = element_text(color = "#5b4f41"),
strip.background = element_rect("#f8f8f8"),
strip.text = element_text(color = "#5b4f41"),
axis.line = element_line(color = "#5b4f41"))
#ggsave(filename = "all_plot_caret.png", plot = last_plot(), width = 12, height = 8, device = "png", dpi = 500) #write an image of the plot
all_pred_pls <- predict(raster_all, fit_all_pls, drop_dimensions = F)
all_pred_feb_pls <- all_pred_pls %>%
filter(date == "2022-02-24") %>%
st_crop(y = boundary, crop = TRUE, as_points = FALSE)
ggplot()+
geom_stars(data = all_pred_feb_pls, aes(x, y, fill = prediction)) +
scale_fill_gradientn(colors = rev(calecopal::cal_palette("desert", 2, "continuous")))+
labs(y = "Northing",#labels
x = "Easting",
fill = "Electrical
Conductivity
(mS/cm)")+
ggtitle("Soil Salinity")+#plot title
theme(plot.title = element_text(color = "#5b4f41", size = 16, hjust = 0.5),#plot themes
plot.background = element_rect("white"),
panel.background = element_rect("#faf7f2"),
panel.grid = element_line(linetype= "longdash", color = "#f0ece1"),
axis.text = element_text(color = "#5b4f41"),
axis.title = element_text(color = "#5b4f41"),
strip.background = element_rect("#f8f8f8"),
strip.text = element_text(color = "#5b4f41"),
axis.line = element_line(color = "#5b4f41"),
legend.title = element_text(color = "#5b4f41"),
legend.text = element_text(color = "#5b4f41"))
mapview::mapview(all_pred_feb_pls)